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Streszczenie 

When interpersonal interactions between individuals are described by the (di- 
screte or continuous) dynamical systems, the interactions are usually assumed to 
be instantaneous: the rates of change of the actual states of the actors at given 
instant of time are assumed to depend on their states at the same time. In reality 
the natural time delay should be included in the corresponding models. 

We investigate a general class of linear models of dyadic interactions with a con- 
stant discrete time delay. We prove that in such models the changes of stability 
of the stationary points from instability to stability or vice versa occur for various 
intervals of the parameters which determine the intensity of interactions. The condi- 
tions guaranteeing arbitrary number (zero, one ore more) of switches are formulated 
and the relevant theorems are proved. A systematic analysis of all generic cases 
is carried out. It is obvious that the dynamics of interactions depend both on the 
strength of reactions of partners on their own states as well as on the partner's state. 
Results presented in this paper suggest that the joint strength of the reactions of 
partners to the partner's state, reflected by the product of the strength of reactions 
of both partners, has greater impact on the dynamics of relationships than the joint 
strength of reactions to their own states. The dynamics is typically much simpler 
when the joint strength of reactions to the partner's state is stronger than for their 
own states. Moreover, we have found that multiple stability switches are possible 
only in the case of such relationships in which one of the partners reacts with delay 
on their own state. Some generalizations to triadic interactions are also presented. 
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1 Introduction 

Interpersonal relationships such as dyadic interactions are very complicated and difficult 
to predict. One of the possibilities of making some predictions is to use mathematical 
modeling. Mathematical modeling of interactions between humans is a relatively new 
branch of applied mathematics. Typically, in mathematical sociology, reactions of indivi- 
duals to stimuli are assumed to be instantaneous, both in discrete and continuous models, 
cf. for example |2l [HI El El [I3l III IIS [111 [111 Ell [22] and the references therein. In 
the basic models of relationships between two partners the linear approach can be used, 
cf. |61 EH It is obvious that such simple models cannot describe the relationships 
precisely. However, they can capture some key features of them. On the other hand, it is 
well known that all natural phenomena are delayed with respect to its stimuli. In the case 
we study, taking into account the delay of reactions can better describe the real dynamics 
of interpersonal relations because of the natural tendency of some people to analyze their 
own and/or partner's reactions and act after that, which necessarily takes some time. 

In the simplest approach, the stimuli are the actual states of both partners. Instanta- 
neous reactions correspond formally to the assumption that the speed, the rates of change 
of the states and the actual states of the partners are considered at the same instant of 
time. However, because time delays are present in human interactions, the reactions to 
stimuli are naturally retarded. Whereas in biological applications the delay has been in- 
troduced many years ago, cf. the Hutchinson model [L2\ which has been proposed in 1948 
as a delayed version of the logistic equation, and the models with delays have been already 
considered in numerous biomathematical literature, cf. for example [H [H [TOl [16] and the 
references therein, the delay in systems modeling human behavior have not received as 
much attention as it possibly should. 

In this paper we study possible scenarios of the evolution of the linear models in which 
the constant in time discrete delay is introduced. We carry out a systematic study of the 
models with delay present in different stimuli terms corresponding to the states of the 
partners. Although in general nonlinear models are more appropriate in the description of 
complicated systems, the mathematical approach we use in the paper allows us to simplify 
the description to the linear one. More precisely, we study small deviations of the partner 
states from some states (which are the equilibrium or stationary states of the system) and 
check whether these deviations tend to zero (stable case) or not (unstable case). Because 
the considered deviations are small, we are able to approximate the nonlinear dynamics 
using linear models, as usually in the stability theory for differential equations without 
delays (ODEs) and with delays (DDEs), cf. e.g. |TOl[TT]. 

In general, different scenarios of influence of a discrete delay on the stationary state 
of the model are possible. 

1. The delay can have no effect on the stability of the stationary state. 

2. Sufficiently large delay can destabilize the stationary state which is stable without 
the delay. 

3. The appropriate values of delay can stabilize the stationary state which is unstable 
without the delay. 
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The second case is the most typical. In this paper we show however, that the third case 
occurs in the systems which are modeled by two linear DDEs for different intervals of the 
coefficients determining the model. 



1.1 Basic Romeo and Juliet model 

As an introductory example, let us consider the model of the mutual interactions of two 
persons, Romeo (R) and Juliet (J), without the delay, cf. |21] [22] . Let the functions 
R{t) and J{t) denote the state of R and J, respectively. The state of the partner can be 
interpreted as the satisfaction, the level of happiness, etc. from their mutual relation. Po- 
sitive values reflect the positive attitude (good emotions) towards the partner, or positive 
experience from the relation, negative values reflect the negative, bad experience. 

In [2T| |22]. the following model has been considered (we cite the author): "R loves J, J 
is a fickle lover. The more R loves her, the more J wants to run away and hide. When R 
gets discouraged and backs off, J begins to find him strangely attractive. R on the other 
hand tends to echo her: he warms up when she loves him, and grows cold when she hates 
him.". In f2T\ this model is described by the system of two linear ODEs 

r m = aj(t), 

1 j{t) = -bRit), ^ ^ 

where a, 6 > 0, according to the interpretation above (other nontrivial cases lead to 
divergent solutions). The solutions to Eqs. ([1]) oscillate around the equilibrium at (0,0) 
which is a center in the phase space (i?, J). 



1.2 General R and J model 

Now, consider more general situation. Let R and J have their fixed goals or preferences, 
ideals of their states in the relation, which we denote by R* and J*, respectively. Assume 
that the rate of change of their actual states depends linearly on two factors: the diffe- 
rences between their own actual states and their own ideals, and the difference between 
the actual states of the partners. Such model has been studied in and reads 

kit) = c[R*-R{t)] + d[J{t)-R{t)l 
Jit) = e[J*-J(t)] + f[Rit)-Jit)]. ^ ' 

The model enables to assign clearly personal features of the persons engaged and com- 
pute the evolution of the system easily. The coefficients c, d, e, / can be positive or 
negative, and describe different psychological traits of the partners, for example contra- 
rian, cooperator etc. We refer the reader to [B] for the discussion of different aspects and 
interpretation of the model. 

We are interested in the temporary evolution of the states of the partners. Let us 
denote by {Req, Jeq) the stationary state of this linear inhomogeneous system. The stability 
properties of this state do not change under the linear transformation of the independent 
variables. Substituting 

r{t) := R{t) - Req, j{t) := J{t) - J,q, 

we obtain the homogeneous system 

r(t) = aiir{t) + ai2jit), 
3 it) = «2ir(t) + a22j(^), 
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with the equihbrium (0, 0), and general coefficients ati G M, A;, / G {1, 2}. 

The conditions for stabihty of the equihbrium for this system can be found in textbooks 
on ODEs, cf. [TT]. In terms of the dyadic interactions they can be interpreted as conditions 
of stable relations. They strongly depend on the signs of the coefficients a^i and the 
relations between them. The choice of signs of the coefficients aki specifies different 
psychological traits of the partners. The case an > 0, ai2 > can be referred to as 
"eager beaver", cf. |22]: "R gets excited by J love for him, and is further spurred on by 
his own affectionate feelings for her." The case an < 0, ai2 > is referred to as "cautious 
lover". The case an > 0, ai2 < could be called "narcissus", and an < 0, ai2 < - the 
doubly cautious, or "stoic lover". It should be noticed that the model is symmetric with 
respect to its variables, therefore the same interpretation as we have given for R is also 
valid for J. 

1.3 Delayed R and J model 

The models of dyadic interactions described by Eqs. ([3]) do not take into account possible 
delays in the reactions of the partners to their instantaneous states. In other words, the 
rate of change of the personal state depends only on the actual state of this person and 
that of the partner. 

In this paper we consider situations in which the delays of the reactions of the partners 
on their own states and/or the states of their partners are taken into account. In general, 
the delay can have various effects on the dynamics of the partnership, for example we 
show that delay can help two cautious lovers to stay in a stable relationship although 
their relationship would lose stability without it. 

We consider the case when the contributions to the rates of changes of the states r{t) 
and/or j{t) can be influenced by the past states of one or two of the partners, that is 



where r^; > 0, A;, / G {1, 2}, are the delays of reactions, that is for r^jt > 0, A; G {1, 2} the 
person reacts with delay on his/her own state while for Tki > 0, k ^ I, k,l ^ {1)2}, the 
delay is present in the reaction to the partner state. It is obvious that if Tki = for some 
k,l G {1,2}, then this type of reaction is instantaneous. 

In general, it is difficult to study the dynamics of Eqs. (jl]) for different values of delays 
T/ci 7^ Tmn foi kl ^ 17171. Howcvcr, the case when the contribution or contributions from 
the past states is delayed by the same time interval can be analyzed almost completely. 
Thus, we consider Eqs. (j3j) with one delay r present in some of the reactions terms, while 
other reactions remain instantaneous. Moreover, it turns out that even in such drastically 
simplified model interesting effects occur, in particular the appropriate delay can stabilize 
the equilibrium outcome of the evolution of the states of the couple. Such effect has been 
also found in [3j for a model with an additional linear term. 

It is well known, cf. [5], that if the delay is large enough, and the equilibrium is 
unstable for at least one value of the delay r > 0, then the trajectories of the considered 
system diverge for delays large enough. In real applications the behavior of the partners 
for bounded delays is more relevant. Let us consider the arbitrary maximal value of delay 
Tmax- We can put the following questions. 

1. The equilibrium is unstable if there are no delays. Can a larger delay < r < Tmax 
stabilize it? 



r(t) 

m 



anr(t - m) + auJit - t^) 
a2ir{t - T2i) + a22j(^ " ^22) 



(4) 
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2. The equilibrium is stable if there are no delays. If the delay is large enough, < 
r < Tinax, it has the stabilizing effect. What conditions on the reactive parameters 
a^; enlarge the interval of stability? 

3. How many stability switches are possible for the situations 1. and 2.? 

4. For given parameters of one partner, what can the second person do in order to 
stabilize the relationship? When is there a possibility to gain stability by changing 
one's attitude towards themselves and the other personSs state, or the relevant 
delays? 

The last question can also be formulated as: can a relationship be in command of one 
flexible person, when the partner is stiff in attitudes and level of refiexivity? 

In the following we address these issues, prove the relevant results and discuss nume- 
rical examples. 

1.4 Notions of stability and instability 

Now, we introduce the notion of stability, instability and stability switches. 

• We say that the equilibrium is stable, when solutions to Eqs. (jl]) tend to this state 
in time. 

• We say that the solution is unstable, when solutions to Eqs. (jl]) diverge. 

• We say that the system is stable when its unique equilibrium is stable. 

• We say that the system is unstable, when the unique equilibrium is unstable. 

• We say that the stability switch occurs if there is the change of stability of equili- 
brium from unstable to stable or vice versa, caused by a change of a parameter of 
the system. In this case typically a Hopf bifurcation occurs. 

For precise definitions we refer to any textbook in ODEs, cf. [llj, or DDEs, cf. |10] . 

It should be noticed that for equations of the form (jlj) there is a unique equilibrium 
(0,0), while for the more complicated, especially non-linear models, there can be more 
equilibriums. In general, we can also consider another type of stability, the Lapunov 
stability, when solutions starting near equilibrium stay near it but not necessarily tend 
to it. Typically, Lapunov stability is non-generic, which means that the small change of 
the model parameters leads to stability or instability of equilibrium in the sense described 
above. In the paper we mainly focus on the generic cases, as typical in analysis of the 
models of natural phenomena. This approach reflects the fact that in reality the models 
parameters do not stay exactly the same all the time but can change in time and therefore, 
non- generic cases are very difficult to observe. 

The terms on the right-hand side (rhs) of Eqs. (jl]) are addressed as stimuli terms or 
reaction terms. 

In the next three sections we consider the models with the delay placed in one, two, 
three and four terms. In section V we introduce a model of triadic interactions. In section 
VI we conclude and discuss some open problems. In Appendix the proofs of theorems 
stated in the paper and some other notes on the topic are presented. 
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2 Single delay 



First we consider the relationships in which only one of the stimuli terms is delayed. There 
are two families of such models. In the first one the evolution of R's state depends on his 
state in the past, that is instead of the general system (jl]) we study 

r(t) = anr{t-T) + aujit), 

j(t) = a2ir(t) + a22j{t). ^ ' 

In other words, R reacts immediately to the state of J, whereas is reflexive to his own 
states, and reacts towards them only after the delay r. This is a reasonable model when it 
comes to a person who tends to analyze the environment rather than themselves and their 
own state affects them only after some deliberation. J is more self-sensitive in this model, 
she tends to update her attitude on the basis of temporary level of private satisfaction. 
Mathematical properties of this model are described by Theorems 1 and 2 below. 

In the second model the evolution of the R's state depends on the J's delayed stimulus, 
whereas his reaction to his own state is immediate, occurs without delay. It is described 
by the following form of Eqs. 

r{t) = anr{t) + aujit-r), 

3 if) = a2ir(t) + a22j{t). ^ ' 

This system describes the situation when R is impulsive about his internal states and 
follows his satisfaction level consciously by reacting to it immediately, but at the same 
time he is cautious about J and reacts to her levels of satisfaction with delay. This may 
result from carefulness as well as lack of knowledge about her feelings. Mathematical 
properties of this model are described by Theorem 3 below. 

It is obvious that exactly the same models can be considered when we replace the 
model variables, that is the roles of R and J are replaced. We omit details for such 
models because its dynamical properties are the same as those presented below. We look 
for stability switches of the equilibrium for both models and We remind, cf. for 
example [3], that if at least one stability switch occur, then (0,0) is unstable for delays 
large enough. Note also that Eqs. ([5]) and (|6]) are stable for r = if 

Q-ii + < and 011022 — 0,12021 > 0. (7) 



Oil O12 

O21 O22 



Inequalities ([7]) can be also expressed by the properties of the matrix A = 
which is formed by the model coefficients, namely 

trA = Oil + 022 < and det A = 011022 ~ O12O21 > 0. 

Note also that for trA = there is Lapunov stability if det A > 0. 
In Appendix we prove the following theorems. 

Theorem 1. For I012O21I < I011O22I 

1. if the system ^ is unstable for t = 0, then it is unstable for all r > 0; 

2. if the system ^ is stable for t = 0, then there is an unique stability switch (for 
which there is the Hopf bifurcation) . 
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Thus, for |ai2a2i| < |fliifl22| there is at most one stabihty switch for all admissible 
reactive parameters aki- 

This result suggests that if the joint strength of reactions of partners to their own 
states, reflected by the product of both strength |aiia22|, is larger than the joint strength 
of reactions to partners' states, there is no stability at all or the system loses stability 
once for all. The conclusion is that, as far as such a relationship was stable and in a result 
of slight moderation of parameters (which usually accompanies personal development of 
partners) starts to shiver, the delayed person should just accelerate reactions and the 
stability will come back to the system. 

Note that the behavior of the system depends continuously on the parameters of the 
model, which implies that decrease in delay will cure the relationship even if it loses 
stability as a result of a slight change in attitude of partners (i.e. some of the parameters 
a^i, k, I G {1, 2}). This situation corresponds to reality and natural development of people 
because it is common to moderate attitudes and reactions in long periods of time. 

Similar interpretation can be given to the results for the model ([6]) , cf . Theorem [3] 
below. 

For 1 0120-21 1 > |fliict22| we define 

A := (oii + 022)^ + 4012021(022 - Oil). (8) 

For A < 0, any of the systems considered in this paper has no stability switches, which 
is proved in Appendix. The case A = is non-generic and therefore we omit it in our 
analysis. For A > and i G {0, 1} we denote 



^ (^al^ - 0^2 - 2012O21 ± Va) , oJi := y/yl for yi > 0, (9) 

(10) 



O12O21O22 

011(^2 _^ q2^) 



Theorem 2. //I012O21I > I011O22I; then for all natural n G N, there exist intervals ofa^i, 
k,l G {1,2}, such that Eqs. ([5]) has n stability switches. The set of sufficient conditions 
reads 

A > 0, al^al2 - 2012O21 > (11) 
and fori e {1,2,. ..,2i\^]-l)} 

1 I (2j-2)7r ^ 1 , (2j-2)7r 

— arccos On + — < — arccos Oi + -^^ —, , , 

1 I (2i-2 TT ^ 1 I (2i)7r \^^) 

— arccos oi + -^^ — < — arccos on + 

The appearance of a stability switch implies instability of the system for r — )■ +00. 
Thus, in the generic case, if the system is unstable for r = 0, the number of switches must 
be even, while if it is stable for r = 0, the number of switches must be odd. However, 
we should remember that in reality the delay is bounded, r < Tmax, and therefore the 
number of switches observed in the interval [0, Tmax] can be different than those obtained 
from the analytical investigation when we let r — +00. 

Hence, for I012O21I > I011O22I, there may be narrow segments of delay values for which 
the relationship is stable. This means the partners cannot know what will be the impact of 
decreasing or increasing the delay on the stability of the relationship, which implies that if 
the system loses stability, it is no longer easy to regain it on the basis of conscious acting. 
This system describes a relationship in which the joint strength of reactions of partners 
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to their own states is weaker than the joint strength of reactions to partners' states. This 
contributes to a common behef that relationship in which partners are more emotionally 
focused on partner's states than on their own, are shaky, because it is essential to hold 
one's own internal source of satisfaction and emotional spine. 

Remark 1. We can also consider other types of linear models, e.g. the situation in which 
the evolution ofr{t) is influenced by its current value r{t) and the past value r{t — T) (and 
the partner's actual state j(t) ) 

( r{t) = aiir(t-r) + ai2j{t) + ai3r{t), 
\jit) = a,,rit) + a,2jit) . 

An example of two stability switches for this model has been considered in JB^. 

We note that for ai3 = the system (1131) is identical to (jSj). Since the stability 
properties of the system depend continuously on the parameters, there exists an interval 
around ais = for which the stability properties of both systems are the same. Since any 
number of switches can be obtained for Eqs. ([5]), the same is true for Eqs. (fT3|) with 
small enough. We argue it more precisely in Appendix. 

For the second model ([6]), in which the rate of change of the R's state depends on the 
J's state in the past, unlike for Eqs. ([S]), at most one stability switch can occur. 

Theorem 3. For |ai2a2i| < |ctiifl22| there are no stability switches for the system 
For |ai2a2i| > |aiia22| 

1. if the system ^ is unstable for r = 0, then it is unstable for all r > 0; 

2. if the system ^ is stable for r = 0, then there is an unique stability switch. 

Remark 2. We note that if we add the delay to the simplest model ([7]j, which is the 
particular case of 

with ai2, 021 > 0, it turns out that, apart from particular initial data, and a measure zero 
set of discrete values of delays, the states of the partners diverge. This system can be 
rewritten as r"{t) = —uj'^r{t — r), where bJ^ = 012021? which only gives stable results for 
T = 2kTi/uj, A; G N. The same is true if both partners react with different delays, Ti and 
T2, respectively, because then the system can be rewritten as r"{t) = —uj'^r{t — (ri + T2)), 
which means it behaves in the same manner, but with the delay r = Ti + T2. However, 
this is stability in the Lapunov sense and therefore it is not preserved for the parameters 
near these critical values. Thus, in this case the delay cannot lead to stabilization of the 
relation. We omit details. 

Below we give an example of the system with a single delay, for which no stability 
switch occurs and the system is stable for every value of delay. Let us consider Eqs. 
with A < and the system is stable for r = 0, which provides a set of conditions 

(a^i + al^Y + 4ai2a2i(a22 - a?i) > 0, 

On + 022 < 0, (15) 
011022 — 012O21 > 0. 
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Conditions ( IT5|) are satisfied e.g. for tlie coefficients an = = 1, 021 = 022 = —2. In this 
example R is an eager beaver, and J is a stoic/double cautious lover, and both of them 
are much more reactive for partner's state than their own. The delay does not affect the 
stability of the system, it affects however the range of time we have to wait for the system 
to stabilize. This means that in such a relationship R can be reflexive to any extent, and 
the relationship is still stable. 

In figures below we present numerical examples of the stability switches and the Hopf 
Bifurcation (HB) for different choices of the parameters of the first model ([5]). In all 
figures the coefficients au of the system ([5]) are represented as the coordinates of the 
vector [ail, ai2,a2i, 0,22] ■ The initial data are r{t) = 1 and j{t) = 1, t G [— r, 0]. 



T=0.35 



T=0.368 



T=0.4 
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Rysunek 1: Time evolution of solutions r(t), j{t) to Eqs. for [—4, 1,-2,— 2] around 
the unique stability switch r ^ 0.3687. Left figure: r = 0.35, stable; middle: r ^ 0.368, 
HB; right: r = 0.4, unstable 



T=0.35 T=0.368 T=0.4 




Rysunek 2: Trajectories of Eqs. (|5]) in the space (r, j) for [—4, 1, —2, —2] around the unique 
stability switch r ^ 0.3687. Left figure: r = 0.35, stable; middle: r ^ 0.368, HB; right: 
T = 0.4, unstable 

In Fig. [T]we show the trajectories r{t), j{t) for delay values close to the unique sta- 
bility switch for the reactive parameters a^/ corresponding to p. 2 of Theorem [H that is 
[011,012,021,022] = [~4, 1, —2, —2]. For delays below the switch value, the equilibrium is 
stable, for larger delays is unstable. The left diagram shows the trajectory which converge 
to the equilibrium. The middle one visualizes the limit cycle appearing as the result of 
Hopf bifurcation, and the right one presents the diverging trajectory. In Fig. [2] we show 
the trajectories in the space (r, j), corresponding to above values of the delay. Simi- 
lar behavior is obtained e.g. for the reactive parameters [—4, —2, —5, —6], [—4,4,4, —9], 
[—5, —2,2, —5]. As it can be seen at least two of o^/, that is On and O22 are negative in 
this case. 

In Fig. [3] we show the trajectories around two stability switches, satisfying the as- 
sumptions of Theorem[2]for the coefficients [5, —4, 3,-1]. For delays values below the first 
switch, the equilibrium is unstable, for delays between the two switches the equilibrium 
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is stable. The left diagrams show the trajectories which converge to the equilibrium. The 
middle one indicate the corresponding limit cycles, and the right ones show the divergent 
trajectories. Similar behavior is obtained e.g. for the reactive parameters [9,4,-6,1], 
[7, 5, —3, 2], [5, —7, 4, 2]. Now, at least two of a^/ are positive. 



i:=0.5 



T=0.522 



1=0.54 





r(t) 

T=0.57 






r-(t) 






r(t) 
T=0.722 







r(t) 

T=0.75 



r(t) 





r(t) 




Rysunek 3: Trajectories of Eqs. for [5, —4, 3, —1] around two stability switches: the 
upper three from unstable to stable equilibrium, the lower three from stable to unstable 
one. Middle diagrams indicate the limit cycles 

In Fig.HJwe show the trajectories around three consecutive stability switches r ^ 0.555, 
T ^ 1.11, T ^ 2.1 for the coefficients [—2, —4, 3, —2]. The left diagram in each row shows 
the trajectory which converge to equilibrium. The middle one approximately visualizes 
the limit cycle, and the right one shows the diverging trajectories. The corresponding 
diagrams which present the time dependence of r[t), j{t) are similar to those in Fig. [H 
and therefore are not shown. Similar behavior can be obtained e.g. for the reactive 
parameters [4, 6, —5, —5], [5, —7, 7, —7], [2, 8, —1, —2]. Here, at least one of a^i is positive. 



3 Delay in two terms 

There are four types of models in which two of the four stimuli terms are delayed. 

r{t) = anr{t-T) + aujit-r), 

j{t) = a2ir{t) + a22jit), ^ ' 

r{t) = aur{t-T) + aujit), 
j{t) = a2ir{t-T) + a22j{t), 

r(t) = anr{t-T) + aujit), 
j{t) = a2ir{t) + a22j{t-T), 

r{t) = anr{t) + auJit-r), 

j{t) = a2ir{t-r) + a22j{t). ^ ' 
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1=0.5 T=0.555 T=0.6 




-10 1 -10 1 -10 1 

r(t) r(t) r(t) 



Rysunek 4: Trajectories of Eqs. for [—2, —4, 3, —2] around three stability switches. 
Middle diagrams indicate the limit cycles 

The first model corresponds to the situation in which R is reflexive towards both himself 
and J. In the second model both partners consider the R's state with the delay. In the 
third model both partners are reflexive towards their own states. Finally the fourth one 
corresponds to the situation in which R and J are reflexive towards their partners. As 
before, we can consider symmetric models with the roles of J and R exchanged, but for 
such models the dynamics is the same. 
We have proved the following 

Theorem 4. 1. Eqs. fT^) . [T?\) and (fT9|) have at most one stability switch. 

2. If ail + (^22 = 0, then Eqs. [W\] have at most one stability switch. 

3. If a22 < 0, (212021 < and an is close enough to 0, then Eqs. f ITS]) can have arbitrary 
number of stability switches. 

For Eqs. f|T8|) with an + 022 7^ 0, it is difficult to analyze stability switches with the 
method presented in Appendix for other systems since the quazi-polynomial has two irre- 
ducible trigonometric terms. However, we are able to provide arguments for the statement 
it does display the same features as Eqs. (|5]) and multiple stability switches are possible. 
We present these arguments in the proof of Theorem H] in Appendix. 
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4 Delay in three and four terms 



There are two families of models in which only one of the stimuli terms is not delayed. 

f(t) = aiir(t-r) + aujit-r), 
j{t) = a2ir{t-T) + a22j{t), 

and 

f(t) = aiir(t-r) + aujit-r), 

3 it) = «2ir(t) + a22j{t-T). ^ ' 

In the first case J reacts immediately only on her own state, in the second case only to 
R's state. 

Finally, there can be the case when both partners react with the delay on both states, 
that is 

r{t) = aiir(t-r) + aujit-r), 
j(t) = a2ir(t-r) + a22j{t-T). 

The quazi-polynomials for systems f l20|) and (12T]) are presented in Appendix. They 
contain irreducible trigonometric terms, as in the case of Eqs. (fTB|) . This means that 
we are not able to formulate appropriate theorems, while for Eqs. (122]) we have proved, 
cf. Appendix, the following 



Theorem 5. // system (l22j) is unstable for r = 0, then it remains unstable for any r > 0. 
// system (122|) is stable for r = 0, then there is an unique switch of stability and the Hopf 
bifurcation occurs. 

The result of Theorem [5] implies that if all reactions of both partners are equally 
delayed, then the dynamics of the system is similar to those for one variable model, 
cf. e.g. [20]. 

Coming back to systems f|20p and f|2ip . we might think analogously to the previous 
models and check the dynamics for possible simplifications. For both models all simpli- 
fications (cf. Appendix) give the same results with at most one stability switch. On the 
other hand, we can provide some arguments that multiple stability switches are possible. 
We present it in Note to systems with three delays in Appendix. However, numerical 
simulations have not bring any other possibilities than no stability switches or only one 
stability switch. We do not present such simulations in the paper because they are almost 
the same as presented for the models with single delay. 



5 Triadic interactions 

One of the possible extensions of the considered models may be obtained by introducing 
the third person. Such generalization has been already suggested in [21], and recently 
investigated in [1]. In particular, in the latter work the authors analyzed periodic and 
quasi periodic solutions for linear and nonlinear models of three-person love relations 
without time delay. 

Let Paris (P) enters the relationship of R and J. There are many possible ways of 
modeling such a triad with delays. We discuss two examples described by the following 
systems of three DDEs 

r(t) = aiir(t) + auj^t), 

jit) = a2ir{t-r) + a22j(^) + a23P(i-^), (23) 
pit) = a^23it) + a33P(i), 
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{r(t) = aur{t) + ai2j{t), 
j{t) = a2ir{t) + a223{t-T) + OasPW, (24) 
P{t) = a32j{t) + a33p{t). 

Both models describe the fact that R and P do not interact directly, therefore do not react 
to the state of the opponent, whereas J needs time to reconsider the situation and react. 
In the model (125]) the pragmatic J analyzes and reacts with delay to the states of both 
men, whereas in the reaction of introvertic J to their states is impulsive, immediate, 
whereas the influence of her states on the speed of change of her actual state is delayed. 

The analysis of the relevant quazi-polynomials which have the third order terms be- 
comes complex. In a particular case of Eqs. (123|) in which R and P react identically 
to their own states: an = 033 = a, we obtain an interesting result. Substituting 
s{t) := a2ir{t) + a23p(t) we achieve 

( j{t) = a22j{t) + l-s(t-r), ^25) 

\ s{t) = (ai2a2i + 023032) j(^) + o-5(^)- 

Thus, we can perceive this situation as if J would interact with an "averaged" man whose 
state is s, and is influenced by this state with the delay. We obtain the model ([6]) with 
the roles of R and J interchanged, in which only one stability switch is possible. Thus, 
when both men have the same attitude towards themselves, there is at most one stability 
switch. 

Below we present an example of numerical results for the second model (125|) for the 
case with three stability switches. Analyzing the relevant quazi-polynomials we can take 
advantage from the fact that stability of the system depends continuously on parameters 
in order to find examples of triads with multiple stability switches. For example, in 
order to obtain a system with n stability switches, we may consider a subsystem with 
only 011,012,021,022 to be non-zero parameters and obtain R-J system with stability 
switches. Then, we may put 032 = 012, 033 = an, 023 = 021 to create a system to be 
symmetrical towards R-P axis. This system will share the number of switches with the 
subsystem. From this point we may try to change values of chosen parameters until 
the number of stability switches jumps to another value. It can be shown that multiple 
stability switches can also occur for adequately chosen intervals of the parameters of the 
models. 

We have used this method to obtain a system with three stability switches for para- 
meters an = —28, ai2 = —74, 021 = 76, 022 = —35, 023 = 76, 032 = —42, 033 = —28. For 
these parameters, the model (12^ describes a situation when after deliberation J reacts 
negatively to her own states, but is quick in response to both R and P, and follows their 
emotions. R and P are similar in reactions and tend to neutralize the relationship: the 
more satisfied J is, the more satisfaction they lose and the opposite. They are the same 
ambivalent towards themselves. The behavior of this system is shown in Fig. [S] 

6 Summary and conclusions 

In the paper we have studied a general family of linear models of dyadic relationships 
with one constant time delay present on the rhs of the systems. We have used the 
mathematical approach presented in [5]. The only linear systems of two DDEs which can 
not be completely analyzed using the presented method are those for which the quazi- 
polynomial has the irreducible term ae~'^'^ + 6e~^'*''^, + 7^ 0, which appears if the delays 
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is present in the 'cross' terms of the rhs of the considered system (since then the auxihary 
function F{y) contains trigonometric terms, and in general analytic considerations are 
very difficult). 

The results presented in the paper, stated mainly in Theorems [T]|U show that the 
joint strength of reactions of partners (where the joint strength means here the product of 
appropriate strengths of both partners) to the partner's state has greater impact on the 
system dynamics than the joint strength of reactions on the partner's own state. This is 
reflected in the fact that for |aiia22| < |cti2fl2i|5 that is for the joint strength of reactions 
of partner's own state weaker than the joint strength of reactions to partner's state, 
the dynamics is typically much simpler. More precisely. Theorem [T] for |aiia22| > |cfi2Ct2i| 
implies that the system ([5]) can have at most one stability switch, while Theorem [2]for the 
same system with opposite inequality can exhibit multiple stability switches. Similarly, 
for Eq. ([6]) Theorem [3] yields no stability switches for |aiia22| > |ai2a2i| and possibility of 
one stability switch for the opposite inequality. 

The fact that stability of the system depends continuously on the parameters leads to 
a general observation that adding a term with the interaction parameter close enough to 
to a given system does not change its stability properties. However, there is a substantial 
difference between systems depending on which term is delayed. Among all the analyzed 
two-person interactions, only those described by the system provide the possibility 
of multiple stability switches. In other words, multiple stability switches are possible 
exclusively when the stimulus of a player to his own state is delayed, and the remaining 



14 



three stimuli are not. Otherwise, i.e. when only the stimulus to the state of the partner 
is delayed, or when at least two stimuli are delayed, then at most one stability switch 
occurs. 

In this paper we have considered the interactions modeled by linear systems. It can 
be of interest to consider nonlinear interactions between the partners. On the other hand, 
one could consider the most general linear model with arbitrary discrete delays. Let a^, 
bk, Cki dk denote the constant reactive coefficients, with for example A; = 1, ...,nrr for 
etc., r^*" > 0, r^-* > 0, r^*" > 0, r^-' > denote the relevant delays, corresponding to 
different interaction terms in the relationship R-J. The general model reads 

-r{t) = a^rit - r^) + ^^^'(^ " ^k)^ 

k=l k=l 

j^j{t) = E ^^^(* - ) + E ^^^'(^ - ^fe')' 

k=l k=l 

and covers all the situations considered in this paper. In particular one could consider the 
situation when the (absolute) values of the reactive coefficients decrease with increasing 
values of the delays, which would mean that the states most remote in the past influence 
less the actual dynamics of the behavior of the partners. 
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7 Appendix 

In the proofs presented below we use the approach developed in [5], cf. also |2D]- If W 
is the quazi-polynomial for the system of DDEs with a single delay r, then W^(A) = 
P(A) + Q{X)e~^'^ , where P and Q are polynomials. If all zeros of W (that is characteristic 
values or eigenvalues of the system) are on the left-hand side of the complex plane, then 
the system is stable. If at least one zero is on the right-hand side of the complex plane, 
then the system is unstable. Because zeros of W depend continuously on the model 
parameters, the change of stability is possible only when some zero crosses imaginary axis 
either from the left to the right or from the right to the left-hand side of the complex 
plane. This is the reason we are looking for purely imaginary zeros of W. Substituting 
X := iu one gets 

W{iuj) = P{ico) + Q{iu)e-''^\ 

Looking for purely imaginary eigenvalues we solve W{iu) = 0, that is P{iuj) = —Q{iuj)e^'^'^'^ 
Taking the norm of both sides and substituting u"^ = y we obtain the auxiliary function 
F defined as 

Fiy)■.= \Pi^^)\'-\Qi^^)\', (26) 

and we see that positive zeros of F give the eigenvalues of the studied system of DDEs. We 
also know that the derivative of F at its zero point determines the direction of movement 
of the corresponding eigenvalues in the complex plane, cf. [SI [SO]- In fact, if this derivative 
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is positive, then the eigenvalue crosses the imaginary axis from the left to the right-hand 
side of the complex plane, while for the negative derivative the movement is in the opposite 
direction. 

Proof of Theorem [1], The quazi-polynomial for Eqs. has the form 

W{X) = - a22A - 012021 + e"^^(aiia22 - Aan). 
The auxiliary function F reads 

Fiy) =y'^ + yial2 - an + 2ai2a2i) + a?2«2i - aii«22- (27) 

If the system is unstable for r = 0, then the corresponding characteristic polynomial, 
which coincides with W{X) for r = 0, has zeros on the right-hand side of the complex 
plane. As r increases, roots of the quazi-polynomial should belong to the left-hand side 
of the complex plane in order to stabilise the system. 

Under the assumption |aiia22| > |ai2a.2i| the function F has only one positive zero yo 
and its derivative 

F'{yo) = 2yo + a^^Ai + 2ai2a2i > 0, 

which means that the transition from the right to the left-hand side is not possible, cf . [20j . 

This implies that if the system is unstable for r = 0, then it remains unstable for all 
r > 0. Moreover, if the eigenvalue enters the right-hand complex half-plane, it remains 
there forever, yielding that only one stability switch is possible. 

If the system is stable for r = 0, i.e. 011022 — 012021 > and On -t- O22 < under the 
assumption of Theorem [H then all characteristic values are on the left-hand complex half- 
plane and can cross to the right-hand complex half-plane which implies destabilisation. 
There exists a critical t^. for which HB occurs and the system is unstable for larger delays. 
More precisely, let cjq = where y^ is the unique positive zero of F . Then the HB 

occurs for the smallest r for which Wiiuj^^ = 0, i.e. for the smallest r such that 



^{W{iujo)) = -Wq - 012O21 - oiia;osin(a;oT) + oii022Cos(a;oT") = 0, 
^{W{iujo)) = -022W0 - ctii«22sin(a;o^) - ctiii^ocos(a;oT) = 0, 

i.e. for Tc determined by the values of sine and cosine 

O12O21O22 . / N COo[l^o + 0-22 ^ '^l'20'2l) 

011(0^2 + ^0) aii(a^2 + ^^o) 

which implies 

^ arccos Oq for On < 0, 



(2J 



27r arccos oq for On > 0. 



□ 



Proof of Theorem [21, If I012O21I > I011O22I, then stability switches can occur only if the 
auxiliary function F defined by Eq. (1271) has two positive roots < yo < yi, i.e. if 



'^22 ~ '^11 + 2012O21 > 0, (022 + 2012O21 - On)^ > 4 ((012021)^ - (011022)^) 
Since F'{y) = 2y + (022 — afi + 2012O21) we have 

F'{yo) = -VA < and F'{y^) = > 0, 
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where A is the discriminant of F, and therefore we deduce that for ujq the eigenvalue can 
move from the right half-plane to the left one, while for Ui - in the opposite direction. 
Eqs. ([2HD yield 



sm Wo-r, 



-a;j{w^+a^2+<^i2'»2i) 



COS[UJjTj) = ttj, 

with Qj = a'^^^u^]^a^ ) J = O5 1- For arbitrary j let > be the first critical delay for 
which Eqs. (129|) are satisfied. 

We obtain two sequences of critical delays for which eigenvalues can cross the imaginary 
axis 

2nT[ , , 

Tjn = Tjo + , n e N, j = 0, 1. (30) 

There are two possibilities 

1. If the system in unstable for r = 0, then at least one zero is on the right-hand side 
of the imaginary axis. Stability switch occurs if the eigenvalue crosses the imaginary 
axis from the right to the left-hand side complex plane through the point iuq. This 
can happen only for the sequence (tjo). Thus, 

Too < no- 

In this case there are at least two switches, because for sufficiently large r the system 
is unstable, cf. [3 ED]- In order to have more switches, the system which is unstable 
for r = should satisfy 

Too < Tio < Toi, 



which gives two corresponding inequalities in (|T2|) . 

2. If the system is stable for r = 0, then it has all eigenvalues in the Ihs of the 
complex plane. If Tiq < Tqq, then at least one stability switch occurs. In order 
to have three switches the following (together with the assumptions listed before) 
conditions should be satisfied 

TlO < Too < Til. 

Analogously one can obtain the set of sufficient conditions for arbitrary number of 

stability switches. 

□ 

Now, consider the possibility of 5 stability switches. Let the system be stable for 
r = 0, then there are 5 switches if 

r A > 0, 

|ctl2Q-2l| > |ctllQ'22|; 

0-11 + 0.22 < A aiia22 — 012021 > 0, 
«ii - O22 - 2012O21 > 0, 

Uq arccos oi < cui arccos Oq < ujq arccos Oi -|- 27ruJo, 
Uq arccos oi + 2nuo < Ui arccos oq + 27ra;i < Uq arccos oi + ^ttuq, 

Five stability switches occur for example for the system ([5]) with the coefficients [—1, 3, —2, 1]. 
The corresponding approximate values of the stability switches are 

Ti ^ 0.92, T2 ^ 1.4, rs ^ 3.3, u ~ 4.2, rg ^ 5.46. 
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The system is stable for r G [0, ri) U (r2, T3) U {t^, t^). 

Note to Remark 1. The quazi-polynomial for Eqs. (fT3|) reads 

W{X) = + (-a22 - ai3)A + (-ctnA + aiia22)e"'^'^ + 013022 - 0-120-21, (31) 

and the corresponding auxihary function F{y) = + (022 + '^13 + 2012021— ctii)y + (fli3a22 ~ 
fli2fl2i)^ — (ctiict22)^- We see that for O13 = both functions are the same as for Eqs. 
This yields that for appropriate values of the model parameters there can be arbitrary 
number of stability switches, as for Eqs. ([5]). 

Proof of Theorem [Si The quazi-polynomial and the auxiliary function F have the form 
W{X) = - (oii + 022)A + 011O22 - oi202ie~^^, 
F{y) =y'^ + (oii + 022)1/ + (011022)^ - (oi202i)^ 

If I012O21I > I011O22I, then -F(O) < and there is one positive zero yo, for which -F'(l/o) = 
> 0, therefore eigenvalues can cross the imaginary axis only from the left to the right. 
Thus, 

1. if the system is unstable for r = 0, then te same holds for all r > 0; 

2. if the system is stable for r = 0, then there exists Tc, for which HB occurs, and for 
T > Tc the system becomes unstable (the critical value be calculated as in 
the proof of Theorem [1]) . 

If I012O21I < I011O22I, then the stability switches are possible only if F{y) has two positive 
roots. The necessary condition afi + O22 < for this situation excludes such possibility. 
□ 

Proof of Theorem 14], 

1. For the models (fT6|) and (fT7|) . the quazi-polynomial is of the same form 

W{X) = A^ - 022A + (-011A + 011O22 - oi202i)Ae"^^, 

and the corresponding auxiliary function F{y) = ?/^ + (o22 — o^i)y — (011022 — 012021)^. 
If 011O22 — 012O21 7^ 0, then the free term of the function F is always negative and 
it has one positive root, which implies there is at most one stability switch. 

2. For the model (fT8|) the quazi-polynomial reads 

W{X) = X^ - 012021 + (-Oil - 022)Ae"^^ + 0110226^^^^ (32) 

If in this particular relationship, partners are focused on their own states from 
the past but with the opposite result, which means On + O22 = 0, the quazi- 
polynomial f l52]) reduces to 

W{X) = A^ - 012O21 + 0110226"^^^ 

which enables the analysis as above. The relevant auxiliary function is of the form 
Fiy) = ~ 2012O21?/ + (012021)^ — (011022)^, and it has two roots for all the para- 
meters, since A(F) = 4(011022)^ is always positive. Hence, 
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a. if 012021 < 011022, the free term of the function F is negative, which imphes the 
function has a single positive root, so that there is no more than a single switch 
possible; 

b. if O12O21 > O11O22, the system never starts from stability, because the conditions 
imply jointly that O11O22 — O12O21 < 0. Thus, if we look for more than one stability 
switch, the condition Too < Tiq < tqi must be satisfied, which implies < ^ < 

where coq = (012O21 — 011022)^, (^1 = (012O21 + 011022)^, but < cjq < cui so that the 
first inequality is never satisfied. 

On the other hand, if On = 0, then quazi-polynomial (152]) reduces to 



W{X) = - O12O21 - 022Ae 



-At 
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which enables the analysis, as before. The auxiliary function F{y) = — (2012O21 + 
^22)?/ + '^i2'^2i has two positive roots under the assumption O22 > 4|oi2||o2i| and 
then multiple stability switches can occur, as in Theorem HI Due to the continuous 
dependence on the parameters multiple switches occur also for on sufficiently close 
to 0. 

3. For the model (fT9|) there is 

W{X) = A^ - (oii + 022)A + 011O22 - e"^'^'^oi202i. 

Here we have F{y) = + {al^ + 022)2/ + (011022)^(012021)^ and respectively A{F) = 
(o^j^oIj)^ + 4(oi202i)^ > for any parameters which means F has two roots. Note 
that if we expect both the roots to be positive, their sum shall be positive, while 
sum of the roots equals — (o^^ + O22) and is always non-positive. Thus, there is no 
more than one stability switch. 

Hence, due to the continuous dependence on parameters the system has a single stability 
switch at most. 

□ 

Proof of Theorem O, The quazi-polynomial has the form 

W{X) = A^ - (oil + 022)Ae"^^ + (011O22 - oi202i)e~2^^ 
and the characteristic equation W{X) = is equivalent to 

^2g2AT _ j^^^^ _^ a22)Ae'^^ + 011O22 - 012O21 = 0. (33) 

Let us denote Ac**"^ =: z, a = — (on + O22) = — trA, b = O11O22 — O12O21 = det A. Then 
Eq. (155]) reads 

z'^ + az + b = 0, (34) 



and therefore 



, — o — y/a^ — 4a6 , —a + a/o^ — 4o6 

Ae^^ = or Ae^^ = . 

2 2 

Looking for purely imaginary eigenvalues one solves 



_ -a - Va^ - 4ab ■ _ -a + Vo^ - Aab 
icuie — — , iu72e — — 
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and obtains two sequences of critical delays at which stability stability switches can occur. 

Now, we need to know the direction of transition of eigenvalues in the complex plane for 
critical values of delays. We use the theorem of implict function to calculate it. Derivating 
Eq. (151|) with respect to r we obtain 



which yields 



dz dz dz 
2z— + a— = ^ — = 0, 
ar dr dr 



dX ^ fdX \ ^ dX X' 
-e^^ + Ae^^ — r + A = 



dr \dT J dr 1 + Ar 

Substituting A = iuk, G {1, 2}, r = Tcr, where Tcr is any of the critical delays, one gets 

d'RX 



dr 



which is positive independently on the values of ujk and Tcr- This means that system [22] 
can have at most one stability switch. 

□ 

Note to systems with delay in three terms. For Eqs. fl20|) the quazi-polynomial has 
the form 

W{X) = X^ - a22X + (-aiiA + aua22)e'^^ - 0120216"^^^ (35) 
and we can reduce it in three ways: 

1. ail = yielding W{X) = A^ — a22A H — 0120216"^'^'^ for which the auxiliary function 
F always has one positive zero; 

2. 012021 = yielding W{X) = A^ — a22A + (— anA + aiia22)e~'^'^ and F has the same 
property as above; 

3. a22 = yielding W{X) = A^ — anAe^'*'^ — ai2a2ie^^'^'^ which has the same form as 
for Eqs. 



In all three simplifications at most one stability switch is possible. 

On the other hand one can rewrite the characteristic equation for A = zo; in the similar 
way as when looking for the auxiliary function F and obtain 

II — — 022^1^11 = llctii^w — aiia22 + 012021 (cos cur — zsinc(;r)|| (36) 

and define the auxiliary function 

Gr{uj) = + (022 ~ '^11)'^^ ^ '^ii'^22 ~ '^i2'^2i + 2aiiai2a2i sinwr + 2aiia22ai2'^2i coscjr 
which zeros give the purely imaginary eigenvalues. We see that 

G{0) = — (aiia22 + '^i2'^2i) < and hm G{ijj) = +00 

implying at least one positive zero of G. We also see that for r = there is Go{uj) = + 
(032 — O'li)^'^ — (aiia22 + 012^21) and Go has one positive zero yielding that for small r > 
only one positive zero exists. However, looking at the graphs of + (032 — ^11)1^^ ~ ^ii'^22 
and a^gi^ii + 2011012021 sin cur — 20110220120^21 cos cur one can easily see that for larger r 
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the number of positive zeros increases and there is arbitrary number of positive zeros for 
sufficiently large r. This shows that there can be more than one stability switch for some 
parameters values. 

For Eqs. (12 ip the quazipolynomial has the form 

W{X) = X^ + (-A(an + 022) - 012021)6"^" + 0110226-2^" (37) 

and we can also reduce it in two ways 

• ^11022 = yielding W{\) = + (— A(aii + 022) — ai2a2i)e^^^ for which F has one 
positive zero; 

• '212O21 = yielding 1^(A) = A^ — A(aii + 022)6^^^ + 0110226"^'^^ and W has the same 
form as for Eqs. (122 p . 

In both cases only one stability switch is possible. However, defining the auxiliary function 
G one can also expect multiple stability switches. 
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